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Abstract 

Using Numerical Stochastic Perturbation Theory within three-dimensional 
pure SU(3) gauge theory, we estimate the last unknown renormalization con- 
stant that is needed for converting the vacuum energy density of this model 
from lattice regularization to the MS scheme. Making use of a previous non- 
perturbative lattice measurement of the plaquette expectation value in three 
dimensions, this allows us to approximate the first non-perturbative coefficient 
that appears in the weak-coupling expansion of hot QCD pressure. 
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1. Introduction 

It is well known that, despite asymptotic freedom, QCD can display non-perturbative 
phenomena even in situations where the system is characterised by a large momentum or 
mass scale Q, Q ^> 1 GeV. For instance, in deep inelastic scattering, power-suppressed 
contributions may arise, of the form ~ Eq^/Q™, where the numerator represents a non- 
perturbative contribution related to a certain quark or gluon condensate of dimension n, 
arising from the Operator Product Expansion, and -Eqcd denotes the typical QCD energy 
scale, of a few hundred MeV. 

A conceptually similar situation arises when Q is replaced by a high temperature T, and 
the observable is replaced by minus the grand canonical free energy density, or the pressure, 
p(T). The Operator Product Expansion gets then replaced by a construction of a low- 
energy effective field theory, which in the case of high temperatures amounts to dimensional 
reduction [T]. The non-perturbative scale -Eqcd gets replaced by that of the effective theory, 
~ ot s T 00], where a s = g 2 /^ is the strong coupling constant. In the effective theory, 
the leading non-perturbative condensate has the dimension n — 3. Therefore, the formal 
weak-coupling expansion of p(T) contains non-perturbative coefficients, starting at 0(g 6 ) 2j. 

While 0(g 6 ) may appear to be an academically high order, recent advances in pertur- 
bative QCD have made its determination an issue of practical importance jlj |3] . Indeed, 
perturbative corrections to the non-interacting Stefan-Boltzmann form of p(T) have been de- 
termined at relative orders 0(g 2 ) C(g 3 ) [7j, 0(g 4 \n(l/g)) 0, 0(g 4 ) 0, 0(g 5 ) PH Ej, 
and 0(g 6 ln(l/g)) |T2j, as a function of the number of colours, N c , the number of massless 
quark flavours, Nf, and, most recently, the chemical potentials, //j, that can be assigned to 
the various quark flavours ^H] (as long as they are small enough compared with the tem- 
perature, \fXi\ <T/g [Tlj). Reaching the next unknown order 0(g 6 ) depends, therefore, on 
the inclusion of the non-perturbative term. Moreover, studies with high-temperature observ- 
ables other than the pressure, in which analogous non-perturbative coefficients arise but at 
lower orders, have shown that their inclusion is in general numerically important, in order to 
reliably determine the properties of hot QCD at physically relevant temperatures [THj 

By now, the first step has already been taken in order to determine the non-perturbative 
0(g 6 ) term: the gluon condensate of dimension n = 3 was measured with lattice Monte 
Carlo techniques within the dimensionally reduced effective field theory in Ref. JTj. Given 
that the other parts of the pressure computation have been formulated in the continuum 
MS scheme, however, this result needs still to be converted to the same regularization ^HJ 
IT§] . The purpose of the present paper is to finalise this task with a certain numerical 
precision. Afterwards, only purely perturbative contributions in the MS scheme remain to 
be determined, in order to complete the expression of p(T) up to 0(g 6 ) (in the notation of 
Ref. [12], Pei remains unknown, while Pe2, Pe3 have recently become available [201 EE])- 

We note that apart from formal interests, the computations outlined above may also 
have phenomenological relevance, in the contexts of cosmology and of heavy ion collision 
experiments. We do not elaborate on these issues in the present paper, however, since the 
current status on these fronts has very recently been summarised elsewhere |21j . 

The plan of this paper is the following. In Sec. 2 we set up the notation and describe the 
overall strategy of the computation. Sec. 3 contains a discussion of our computational tool, 
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Numerical Stochastic Perturbation Theory in a covariant gauge. Sec. 4 contains the data 
analysis and our results, while Sec. 5 draws some conclusions. 



2. Basic definitions and overall strategy 



Let us start by considering three-dimensional (3d) pure SU(iV c ) gauge theory in dimen- 
sional regularization. The Euclidean continuum action can be written as 

1 



S 
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d a xC E , C E = — j Tr [F kl F kl \ 

Z 93 
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where d — 3 — 2e, g\ is the (dimensionful) gauge coupling, k, I = l,...,d, F ki = i[D k ,Di\, 
D k = d k + iA k , A k = A B T B , T B are the Hermitean generators of SU(iV c ), normalised as 
Tr [T A T B ] = 5 AB /2, and repeated indices are assumed to be summed over. Leaving out for 
brevity gauge fixing and Faddeev-Popov terms, the "vacuum energy density" is defined as 

1 
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J VA k exp^-^E 
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where V is the <i-dimensional volume, DA k a (gauge-invariant) functional integration mea- 
sure, and we have assumed the use of the MS dimensional regularization scheme to remove 
any 1/e poles from the expression. 

The physical significance of f^jg for high-temperature QCD enters through the framework 
of dimensional reduction [UI22]- Indeed, after the replacement g\ — > g 2 T[l + 0(g 2 )}, 
appears in the QCD pressure as an additive contribution, Sp(T) = — Tf^g |2J [TT]. The 
dependence on regularization (in particular, on the MS scheme scale parameter Ji) disappears, 
once all other contributions of the same order jT2] have been added. 

Now, dimensional reasons and a perturbative computation of ultraviolet divergences 
show that the structure of f j^g is (after letting e — > 0) 
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where d A = A^ c 2 — 1. The non-perturbative constant Bq, which actually is a function of 
N c , is what we would like to estimate in the following. For future reference, we note that a 
logarithmic derivative of f j^g with respect to g 2 immediately produces the gluon condensate: 
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We now go to lattice regularization. In standard Wilson discretization, the lattice action 
5*l, corresponding to Eq. (JTJ), reads 



Sl 



(5) 



k<l 



where n w (x) = Re[Tr E/ju(x)]/JV c , U k i(x) = C/ fc (x)^(x+ Aj)E/J(x + l)U}(x) is the plaquette, 
U k (x) is a link matrix, x + A; = x + ae^, where a is the lattice spacing and e k is a unit vector, 
and 
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Note that the gauge coupling does not get renormalised in 3d, and the parameters g% appear- 
ing in Eqs. (JTJ), (JHJ) can hence be assumed finite and equivalent. The observable we consider 
is still the vacuum energy density, Eq. (J2J), which in lattice regularization reads 



f L = - Jim - In 



(7) 



where T>Uk denotes integration over link matrices with the gauge-invariant Haar measure. 

Being in principle physical quantities, the values of f^jg and must agree, provided that 
suitable vacuum counterterms are added to the theory. Due to super-renormalizability, there 
can be such counterterms up to 4- loop level only [21], and correspondingly 

Af EE f L -f m 

= C 1 -f\n^ + C[)+C 2 ^ + C 3 ^ + C,g!(\n— + C' 4 )+0(gla), (8) 
a \ a 9z J a a \ a A* / 

where the Cj are dimensionless functions of N c . The values of Ci, C 2 , C3, C 4 are known, as 
we will recall presently; C[ is related to the precise normalisation of the Haar integration 
measure and has no physical significance; and C 4 will be estimated below. 

Correspondingly, the gluon condensates, i.e. the logarithmic derivatives of fj^g, with 
respect to gf , can also be related by a perturbative 4- loop computation. Noting that three- 
dimensional rotational and translational symmetries allow us to write 



(10) 



and employing Eqs. (J3J), (JBJ), leads to [T7] 

8 H B « - ^A {i - n ^ - [t + w + w + c 0] } ■ 

The values of the constants trivially related to those of Ci,...,C 4 in Eq. (JBJ): 

d = Ci/3, c 2 = -2N c C 2 /3, c 3 = -8A^ c 2 C 3 /3, c 4 = -8A^ C 3 C 4 and, in particular, 

c 4 = C 4 -i -21n(2JV c ) . (11) 
For N c = 3, the constants read jT7 ] 125 ] I2T) ] 177] 

ci = ^ « 2.66666667 , (12) 

c 2 = 1.951315(2) , (13) 
c 3 = 6.8612(2) , (14) 

Ci = 2.92942132. (15) 



(4tt) 4 V12 768 
Moreover, lattice measurements J7J have shown that, for N, 
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In order to extract Bq, as is our goal, we need to determine the unknown constant c 4 in 
Eq. (fTHJ). This can be achieved by repeating the same setup as above, but by regulating 
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infrared (IR) divergences through a mass regulator, m, instead of confinement. Indeed, the 
difference in Eq. (JHJ) is IR insensitive, and does not change. Thus we can extract C' A this way 
and, from Eq. (fTTj) . c' 4 . We denote quantities computed with a mass regulator with a tilde. 
With a mass regulator, the continuum computation produces 
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where the lower order terms omitted vanish for m — > 0, and B^ depends on the gauge 
parameter a, since the introduction of a mass breaks gauge invariance. If, on the other 
hand, we carry out the same computation in lattice regularization, we expect 
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where the lower order terms omitted go over to the ones in Eq. (jHJ) for m — > 0. Given these 
forms and the IR insensitivity of the difference in Eq. (JHJ), we get 
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where the lower order terms omitted agree with Eq. (jHJ. Comparing the C(^f) terms with 
Eq. (jHJ), we can read off C' A . Subsequently, Eq. (fTTj) gives c? 4 . Inserting finally the result from 
Eq. ()16)) . we arrive at the master relation 



B G = 10.7 ± 0.4 - B L (a) + B m (a) + 
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It remains to determine Bj^(a) and Bi(a). Since the result of Eq. (j2~U|l is gauge indepen- 
dent, we choose the covariant Feynman gauge (a = 1) in the following. A 4- loop continuum 
computation, described in Ref. leads to a small set of fully massive master integrals 
that are known with high precision [21], an d produces (see also Ref. |29|) 



-2.16562591949800919016... . 



(21) 



On the other hand, taking a logarithmic derivative with respect to g\ from Eq. (|18j) . we 
obtain, in complete analogy with Eq. (|10jl . 



up to 4-loop 
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Our task in the following is to compute the right-hand side of this equation for a = 1 and, 
afterwards, to insert the result into Eq. (j2Hj) . in order to estimate Bq. 
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3. Numerical Stochastic Perturbation Theory 



To carry out the limit in Eq. (}2*2*|) requires a 4-loop computation in lattice perturbation 
theory. As the master integrals that appear in higher loop computations in lattice regu- 
larization need to be evaluated numerically in any case, we choose to carry out the whole 
computation numerically This can be achieved through the use of Numerical Stochastic 
Perturbation Theory (NSPT) 30 , pioneered in recent years by the Parma group; a full 
account of the method can be found in Ref. [3~T] . 

In its "purest" form, NSPT can be applied without either gauge fixing or masses as 
IR regulators. Since we compare with a dimensionally regularized gauge fixed continuum 
computation with a mass as an IR regulator, however, we need to introduce the same tools in 
NSPT. The first two subsections describe our general implementation, and the third collects 
some technical details of the computation. 

3.1. NSPT in a covariant gauge. NSPT relies on Stochastic Quantization [32] (for an ex- 
tensive review, see Ref. [33J). In this approach quantum fields are given an extra coordinate, 
r, which is to be regarded as a stochastic time in which an evolution takes place according 
to the Langevin equation. This is in close analogy with the "time" evolution of the Markov 
chain that is used in standard Monte Carlo simulations; indeed, Stochastic Quantization can 
also be used for Monte Carlo simulations [SI] (for a concise review, see Ref. |35j ) . 
For lattice gauge theories, the Langevin equation reads 



where we assume the use of lattice units (a = 1) in the spatial directions. The derivative 



V fcjX is defined [HI as V fc , x = T A V^ X , with V£ >x S[U k (x)] = lim e _ >0 {5[exp(zeT A )f/ fc (x)] - 



5[f/ / fc(x)]}/e, where T A are the generators in the fundamental representation, normalised 
as before. Moreover, r)k is a gaussian noise in the adjoint representation, r/k = r/ A T A . In 
Eq. ([2*3]) we adhere to a precise notation in which the dependence of the solution on the 
stochastic noise is explicitely shown. Since this notation is a bit pedantic, we will drop 
it in the following. It is worth stressing that the evolution dictated by Eq. (|2*3*|) preserves 
unitarity, and that Vfc jX is consistent with partial integration over the Haar measure. 

The main assertion of Stochastic Quantization is that the path integral correlation func- 
tions of the field theory, computed with the Haar measure, can be traded for stochastic time 
averages in the asymptotic r — > oo limit: 



where O is some observable. NSPT is obtained by expressing the solution of the Langevin 
equation as a power series in the coupling constant and by numerically integrating the 
hierarchy of equations that results from inserting this expansion into Eq. (|2*3*|) . In our notation 
the expansion reads [3U1 OH] 




(23) 
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in which N is the highest order one wants to reach in the computation. In our case we need 
to expand the field up to /3~ 4 , that is iV = 8. Note that since Eq. (J23|l guarantees unitarity 
by construction, the terms (x) will automatically inherit the corresponding properties. 

While the spirit of Parisi's and Wu's original paper [32] was to offer a possibility for 
performing perturbative computations without gauge fixing, gauge fixing can naturally be 
added to the framework through the Faddeev-Popov mechanism, like to lattice gauge theory 
in general. The partition function is written as 

Z = / Wl de t A PP exp(- 5l - SoF ), (26) 
where Sl is the Wilson action in Eq. (JHJ). For Sqf, we choose the standard covariant form, 

P 
AN c a 



= J-£fe^(x)f, (27) 



where a is the gauge parameter; d% is the left difference operator; we have defined Lie 
algebra valued fields <pk through Uk = exp(i0fc); and we have written (pk = 4> A T A - The naive 
continuum limit is obtained through the identification 0^ = ag^A^. Once we enforce the 
expansion in Eq. (|25|) for the links, the Lie algebra valued fields are expanded as well: one 
simply needs the Taylor series for 0^ = — iln[l + SUk). 

Note that in standard lattice perturbation theory Eq. (J26j) is further modified by writing 
the functional integration in terms of On the contrary, we directly expand in terms of 
the link variables Uj. (cf. Eq. ()25|l). in terms of which Eq. (|23|) is formulated. Therefore we 
do not need to add any "measure term" to the action. Nevertheless, the values of various 
operator expectation values remain identical, order by order in to those in the standard 
perturbative framework. 

The Faddeev-Popov operator corresponding to the gauge function in Eq. (|2*7|) reads App = 
—d^D). [<j>]. Since Sqf is expressed in terms of Lie algebraic fields, Dk[4>] comes from evalu- 
ating the response of the field <pk to a gauge transformation. The latter is defined in terms of 
the link variables Uj., so it does not come as a surprise that Dk[4>] is not expressed in closed 
form. It can however be written as a perturbative expansion, whose first terms read (see e.g. 
Ref. j2Z], whose notations we follow) 

4M(x) = [i + •. tW _^ w _^ w _^ !w+0 (^)" 



9f+»*(x), (28) 



where dj} is the right difference operator, while is the field in the adjoint representation, 

$ fe (x) = <^(x) F A , (F A ) BC = -if ABC . (29) 

In order to set up the proper Langevin equation, we finally rephrase the Faddeev-Popov 
determinant as a new contribution to the action, det App = exp(Tr In App). 

3.2. Mass regulator. We still need to add a mass regulator to the gauge-fixed framework. 
We note that the same mass is given to the gluon and the ghost fields in the continuum 
computation we want to match to. We therefore modify the total action to become 

S = Sl + Sqp + Spp . (30) 
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The gluonic action has been modified by a mass term, 

Sl = S l + ^^£^(x)^(x) , (31) 

X 

and the Faddeev-Popov action reads 

S FP = -Tr lnA F p , A FP = -d^D k [(f)} + (am) 2 . (32) 

Even if we never make use of ghost fields in our approach, Eq. (}3*2*|) amounts to giving them 
a mass, as is clear from writing the operator App in Fourier-representation. 

3.3. Some technical implementation issues. To treat the Faddeev-Popov determinant 
as a part of the action means that, because of the Langevin equation, one has to face the 
quantity V k}X S FP = -V fc , x Tr [In A FP ] = -T A Ti [VjJ x A FP App], with App from Eq. (JH2J). 
We follow the procedure proposed in Ref. 38] (see also Ref. [34]). Our study is actually the 
first practical implementation, but the essential ingredients are the same as for the treatment 
of the fermionic determinant in unquenched NSPT, as discussed in Ref. |31j . 

By introducing an extra gaussian noise £ normalized as {£,l£, n )(, = 5km one can substitute 1 

-Tr [V£ x App A F p] — > -(Re{a(V^ x A F p) M (App)^ n }) 5 . (33) 

Taking the real part would be unnecessary after the average (...)(, but we find it convenient 
to impose it already before the averaging. The advantage of the form introduced is made 
clear by rewriting 

£*(V£ x A F pWApp)^ n = £*(V£ X A FP )^ , (34) 
where ip = App£. In NSPT, then, we need to compute 

V>W = (Ap p )« e • (35) 

It is worth stressing that the noise £ has no power expansion, while the field ip (like any other 
field in NSPT) is expanded because of the power expansion of App (which is a function of 
the fields 0, i.e. of the fields U). 

That Eq. (J3*5j) can be evaluated efficiently within the NSPT framework stems from the 
expansion of the operator A F p. Once a generic matrix M is given as a power expansion, 

oo 

M = M (0) + P'^M^ , (36) 
the expansion for its inverse reads 

oo 

M- 1 = [M^f' + J^^HM-'f ■ (37) 

i=l 



Here k,l,n should be regarded as multi-indices, including space coordinates and colour. 
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Here the non-trivial terms are obtained through a simple recursive relation: 
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In our case this leads to 
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-i 
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(39) 



j'=o 



Eq. ()39|) states that there is no actual matrix inversion to take. Indeed, [Ap p ] is indepen- 
dent of the <fi fields and its expression is well known: it is the (would-be) ghost propagator, 
diagonal in Fourier space. Note that the mass regulator makes it well-defined at every value 
of the momentum. As for the various orders they are naturally computed by iteration. 
At every order only one application of [App] is needed: this propagator operates on a sum 

of already computed quantities (the lower order ip's). While [App] is diagonal in momen- 
tum space, all the other operators are almost diagonal in configuration space. This suggests 
the strategy of going back and forth from Fourier space via a Fast Fourier Transform. It 
remains to point out that also the expression for x App (and its power expansion) is sub- 
stantially local, so that the big inner product in Eq. (|34j) is not too difficult to deal with. 
Finally, Spp and App are naturally written in the adjoint representation, so that one has to 
devise an efficient way of dealing with cascades of commutators of the fields. 

In order to solve Eq. ()23|) numerically, the stochastic time variable r needs to be discretised: 
t = na T , where n is an integer. We use different values of a T , average over each thermalised 
signal, and then extrapolate in order to get the value of the desired observable at a T = 0. 
Eq. (j23J) is discretized in the standard way jSl] which automatically preserves the unitarity 
of our degrees of freedom: 



where 



f4(x, (n + IX) = e~ lF ^ na ^ U k (x,na T ) 



F k (x,na T )[U,r)] = a T V k;X S[U] + vX%(x,na r ) 



(40) 



(41) 



2(5 X yO mn . 



and we have assumed the normalization {r)(x,ma T )r](y,na T )) v 

We note that Eq. (J4*U|) is only accurate to first order in a T . As a consequence, if the action 
is written as a sum (S = Sl + Sgf + Spp), one can to the same accuracy rewrite Eq. (j4*Uj) 



am 


exact 


NSPT 


0.10 


2.6579 


2.6581(6) 


0.15 


2.6481 


2.6476(6) 


0.20 


2.6379 


2.6379(8) 


0.40 


2.5713 


2.5710(8) 


0.60 


2.4714 


2.4710(8) 


0.80 


2.3485 


2.3481(7) 


1.00 


2.2119 


2.2117(6) 


1.20 


2.0689 


2.0681(5) 


1.40 


1.9251 


1.9247(5) 


1.60 


1.7848 


1.7846(4) 



Table 1. Comparison between "exact" and NSPT results for the coefficient 
of j3~ l for the observable in Eq. ()45|) . The lattice extent is L = 8a, except for 
am = 0.10, where it is L — 10a. 

as 

f4(x, (n + l)a T ) = e"^^ { e - iF P™ U k {^na T ) } . (42) 

The advantage of this form is that from the point of view of program implementation, it is 
easier to first evolve the field by 

FM[U,r)) = a T V k , x S L + v^7% (43) 

(i.e. the contribution coming from standard Wilson action plus gaussian noise) and then by 

F^[U}= a T V k , x (S - S L ) . (44) 

Indeed, F' 2 1[£7] depends on the C/'s only through the 0's, i.e. the first step is performed in 
terms of the U fields, the second in terms of the fields. It is also easy to realize that the 
first step can be implemented as a sequential sweep through the lattice, while the second 
one requires the construction of a global contribution (the inner product in Eq. (J23J). 

4. Data analysis and results 

Since this work is the first time that the Faddeev-Popov procedure was implemented 
in NSPT, much attention was devoted to reliability checks, which we describe in the first 
subsection. The second subsection is devoted to carrying out the limit in Eq. (|2^jl . All the 
numerical values shown in the following were obtained with iV c = 3. 

4.1. Consistency checks. The first checks were performed against the theory without 
gauge fixing, both in 3 and in 4 dimensions. In other words, we set m = 0, and checked that 
we reproduce gauge invariant results for the plaquette expectation value, irrespective of the 
gauge parameter a used. In particular, we made sure that, for fixed volumes (V = L 3 ), we 
could reproduce the 3d results in Ref. [26 , up to 4- loop order. 
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Figure 1. Examples of finite- volume values (open diamonds) and infinite- 
volume extrapolations (open squares) for the coefficient of /3~ 3 (left) and f3~ A 
(right) in the mass regularized plaquette expectation value, Eq. (|43j) . 

We then plugged the mass terms in. A first test was that the leading 0{f3~ 1 ) contribution 
to the plaquette expectation value appearing in Eq. ((22} , 



fr \ _ 1 d A 



1 - 



f \2 N 

,am) \ ^ 



AT3 



[EUi 4sin \-nn,/N) + {amf\ 



(45) 



where N = L/a, was correctly reproduced by the NSPT numerics. Tabled shows the results 
for the coefficient of and gives an idea of the size of errors. It also lists all the mass 
values that we will use in the following. 

As a final preparation we checked, as already explained in Ref. jSH], that if the volume is 
kept finite, and a fit in (am) 2 is performed for measurements carried out with the masses 



shown in Table ^ then the intercept with the axis am 
the massless theory, obtained without gauge fixing j2H|- 



agrees with previous results in 



4.2. Detailed analysis. After the consistency checks, we turn to the actual analysis of 
Eq. (122)1 . It consists of two steps. First, for any given value of m, we need to carry out an 
infinite-volume extrapolation. Second, the infinite-volume extrapolations need in turn to be 
extrapolated to m — > 0, as dictated by Eq. (j221)- Both of these extrapolations are rather 
delicate so let us describe the procedure that we adopt in some detail. 

For the infinite-volume extrapolations, we expect that the dependence on L is exponen- 
tially small, once the volume is large enough. Concretely, inspecting the 1-loop expression 
in Eq. ()45|) for am 1 but at finite volumes, suggests that for mL ^> 1 the behaviour is 

<l-Ili 2 ) L ~7o + ^ 
have all the same sign. However, for our larger masses am > 0.40 (cf. Ta- 



7i 



exp(— mL) + 7 2 exp(— v^mL) + 



(46) 



where 71,72, 

ble[TJ), the volume dependence appears in fact to be dominated by discretization effects not 
contained in Eq. (j4*U|) . Moreover, at higher loop orders other structures also appear and it is 
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1 



0.8 1.0 
1/mL 

Figure 2. An example of different infinite- volume extrapolations for the 
coefficient of /3~ 4 in the mass regularized plaquette expectation value, for 
am = O.fO. The open square denotes the extrapolation described in the text; 
the closed circle is an extrapolation according to Eq. (|%5Jl (with unconstrained 
coefficients) ; and the closed triangle assumes yet another phenomenological fit 
form, So + exp(— mL)[5i + 82 (mi) -4 ]. 

not clear a priori how large mL has to be for them to remain negligible (note that for the 
important small masses am < 0.20 we are only able to go up to mL ~ 2, cf. Fig. |TJ). 

For these reasons, we adopt a practical procedure in the following whereby we increase 
the volume until no volume dependence is seen within the error bars, and then fit a constant 
to data in this range. To be conservative, the resulting error bars are multiplied by a 
factor two. The original "raw" data at finite volumes, and the corresponding infinite-volume 
extrapolations obtained with the recipe just described, are illustrated in Fig. [T] for a few 
masses. However, we have also tried other procedures, like a fit according to Eq. or to 
different phenomenological forms ruled by decaying exponentials multiplied by polynomial 
prefactors. Fig. 121 gives an idea of the effects of these variations on the infinite-volume 
extrapolations. It can be seen that our doubled error bars can cover all the variations. 

Given the infinite-volume extrapolations, we can carry out the extrapolation am — > 0. 
Motivated again by a 1-loop analytic computation (Eq. ()45j) for iV — > 00), we use an ansatz 
allowing for any positive powers of am. There is the problem, however, that the data points 
are more precise at larger masses: the absolute errors decrease roughly as ~ l/(am) 2 , i.e. 
vary by two orders of magnitude. Thus large masses tend to dominate the fit, while the most 
important region should be that of small masses. 

We confront this situation in the following way. First of all, we allow for a high-order fit 
polynomial, and monitor the stability of the results, and the x 2 - y alue of the fit, with respect 
to the order of polynomial, as well as the number of masses that are taken into account. 
Second, we consider both the regular x 2 -function, and modified ones where the errors are 
weighted by am or by (am) 2 , so as to make the contributions of the different data points more 
balanced. All of these fits generally show an extremum as the order of the fit polynomial is 
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Figure 3. The am —>■ extrapolations (dashed lines, with the intercepts 
with the axis am = shown with closed squares) through the infinite-volume 
extrapolated data points (open squares), for the coefficient of [3~ 3 (left) and 
j3~ A (right) in the mass regularized plaquette expectation value. 



increased, with the extremal values differing by less than the statistical errors. Examples of 
the fits are shown in Fig. El 

In order to benchmark this strategy, let us first apply it at 1-loop, 2-loop, and 3-loop 
levels. We obtain 



lim < lim (l — 

rrwO [L->oo x ' L 

lim < lim (l — IIit^t 
lim < lim (l — I1i2)t 

m^0 L->oo x ' L 



2.672 ±0.008 
1.955 ±0.016 
6.83 ±0.10 . 



(47) 
(48) 
(49) 



These numbers are to be compared with the known results in Eqs. fll2j) - (fT4*j) : we find perfect 
agreement within error bars. 

We then repeat the same procedure at 4- loop level. As shown by Eq. (J22J), the extrapolation 
am — > can only be carried out after the subtraction of the logarithmic IR term. The fit 
shown in Fig. El produces 

,2\ 2 



Bl(1) 



2^ 
~27 



x (25.8 ±0.8) = 13.8 ±0.4 . 



Inserting into Eq. (|2~Uj) . a significant cancellation takes place, and we obtain 

5 G = -0.2±0.4( MC )±0.4( NSPT ), 



(50) 



(51) 



where the error "MC" originates from the lattice Monte Carlo simulations in Ref. and 
the error "NSPT" from the analysis in the present paper. Eq. (joT|) is our final result. Let 
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us also record the N c = 3 values for the constant C 4 in Eq. (JHJ) and c' 4 in Eq. (fTUj) . 

C 4 = 10.9 ± 0.3( NSPT ) , c' 4 = 7.0 ± 0.3( NSPT ) . (52) 

5. Conclusions and perspectives 

We have demonstrated in this paper the feasibility of determining the non-perturbative 
constant Bq, defined through Eq. @, by combining previous lattice Monte Carlo results [T7] 
with a 4-loop perturbative matching step. The matching involves a comparison of a contin- 
uum MS computation j2S] with the corresponding lattice regularized computation. The latter 
we have carried out with the help of Numerical Stochastic Perturbation Theory (NSPT). The 
final estimate for the new matching coefficient (with two different conventions) is shown in 
Eq. (jo^j) . Taking into account the Monte Carlo results, the final estimate for Bq is shown 
in Eq. (loTl) . 

We note that within the current errors, Bq is consistent with zero. This is a matter of 
conventions, however; for instance, had we not made the arbitrary choice of including the 
factor 2 inside the logarithm in Eq. ©, the corresponding constant would be non-zero by a 
significant amount. 

Given that the physical pressure of hot QCD is numerically fairly sensitive to Bq [12], it 
would of course be desirable to improve on the accuracy of Bq, both on the lattice Monte 
Carlo and on the NSPT sides. For instance, it would be interesting to repeat the current 
study with traditional techniques [2Zj. Moreover, it should in principle be possible to carry 
out the matching leading to Eqs. (|52*|) by using a finite volume rather than a mass as an 
infrared regulator; for this approach the NSPT side exists already [21], but the 4-loop MS 
computations of Ref. [22] would have to be repeated with techniques discussed at 1-loop 
level for instance in Refs. [4*U] . 

Apart from these challenges, there is now an ever more compelling case for determining the 
last remaining purely perturbative 0(g 6 ) term, i.e. the constant denoted by /?ei i n Ref- [J2|- 
Only after this has been added does the //-dependence of Eq. (jSJ) get cancelled, such that 
the full physical pressure is scale- independent, as it has to be. 
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